Molecular conduction: paradigms and possibilities 
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We discuss the factors that determine the overall shape and magnitude of the current- voltage (I-V) 
characteristics of a variety of molecular conductors sandwiched between two metallic contacts. We 
analyze the individual influences of the contact geometry, the molecular chemistry, the electrostatics 
of the environment, and charging on molecular conduction. The theoretical predictions depend 
sensitively on the experimental geometry, as well as on the theoretical model for the molecule and 
the contacts. Computing molecular I-V characteristics will thus require theoretical understanding 
on several fronts, in particular, in the scheme for calculating the molecular energy levels, as well as 
on the position of the contact Fermi energy relative to those levels. 



I. MOLECULAR CONDUCTION: WHAT IS THE 
UNDERLYING PHYSICS? 

Recently several researchers have measured charge 
transport in single or small groups of organic molecules 
connected to metal contacts [1-9]. In parallel, there 
have been theoretical attempts at understanding molecu- 
lar conduction, both at the semi-empirical [3,10-16] and 
first-principles [17-23] levels. Understanding molecular 
conduction is challenging, since it involves not just the 
intrinsic chemistry of the molecule, but extrinsic factors 
as well, such as the metal-molecule bonding geometry, 
contact surface microstructure and the electrostatics of 
the environment. The aim of this article is to discuss the 
various physical factors that influence molecular current- 
voltage (I-V) characteristics, and our attempts to model 
them both qualitatively and quantitatively. 

A typical two-terminal molecular I-V looks like Fig. 1, 
often with a clear conductance gap [1]. How does one 
understand such an I-V? The first step is to draw an 
energy-level diagram, as in Fig. 2. An isolated molecule 
has a discrete set of energy levels, with a highest occu- 
pied (HOMO) and a lowest unoccupied molecular orbital 
(LUMO), separated by a HOMO-LUMO gap (HLG). On 
connecting the molecule to metallic contacts, two changes 
happen: (i) the discrete molecular levels broaden into a 
quasicontinuum density of states (DOS) due to hybridiza- 
tion with the metal wave functions. Often the DOS re- 
tains a distinct peak structure, in which case it is still 
useful to think in terms of broadened molecular energy 
"levels" ; (ii) the difference in work functions between the 
molecule and the metal causes charge transfer and band 
alignment between the two materials. The molecule equi- 
librates with the contact with an overall chemical poten- 
tial set by the metal Fermi energy Ep, typically lying 
inside the HLG. Under an applied bias the molecule tries 
to equilibrate simultaneously with both contacts with 
bias-separated chemical potentials /ii,2, and is thereby 
driven strongly out of equilibrium. As long as the bias 
is small and /zi,2 he in the HLG, the HOMO levels stay 
filled and LUMOs are empty and there is no current. 



However, when the bias is large enough that either /ii 
or [i2 crosses a molecular level -Emol, that level is filled 
(reduced) by one contact and emptied (oxidized) by the 
other (Fig. 2a) and therefore starts conducting current 
[24]. For opposite bias, the same level starts conduct- 
ing when crossed by the other contact chemical poten- 
tial (Fig. 2b). The net result is that for a spatially 
symmetric molecule with symmetrically coupled contacts 
the total conductance gap is given by ~ 4 (Ep — Emol)- 
Molecular conduction thus depends on both the intrinsic 
molecular chemistry through Emol and the contact mi- 
crostructure through Ep. 
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FIG. 1. Generic molecular I-V and parameters controlling 
it. The current rises when a molecular level -Emol is crossed 
at a bias V m 2(Ef — Emol), where Ef is the contact 
Fermi energy. The overall current magnitude is controlled 
by r c ff = riLV/pTi + IV], Fs being the broadenings of the 
molecular levels by hybridization with the contacts. The cur- 
rent rises over a voltage width set by the thermal broadening 
ksT and by Y = Ti + IV The current is dragged out further 
by the presence of a Coulomb charging energy Uo- 

The intrinsic chemistry of an isolated molecule can be 
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handled with sophisticated quantum chemical codes that 
can be purchased or even downloaded from the Internet. 
Given an appropriate basis-set (for instance, a minimal 
STO-3G basis) and an appropriate model for electron- 
electron interactions (based on first-principles density 
functional, Hartree-Fock or semi-empirical Huckel meth- 
ods), such a code starts with an initial guess density ma- 
trix p to obtain a Fock matrix F (Fig. 3a). It then 
fills up the corresponding energy eigenstates with a given 
number N of electrons according to equilibrium statis- 
tical mechanics to reevaluate p, recalculates F and so 
on, until self-consistent convergence. Our molecular sys- 
tem differs in two ways: (i) it is open, with a vary- 
ing, fractional occupancy of electronic levels; (ii) it is 
trying to equilibrate under bias with two different con- 
tact chemical potentials and is therefore driven strongly 
out of equilibrium. To solve this problem, we mod- 
ify the above self-consistent scheme, as shown in Fig. 
3b. The initial step, solving for the Fock matrix F, is 
kept unchanged from the usual prescriptions in molecu- 
lar chemistry. In this step, one can use semi-empirical 
tight-binding/Huckel-based methods [25], or exploit the 
sophisticated numerical prowess of a standard quantum 
chemical package such as GAUSSIAN'98 [26] to employ a 
density functional theory (DFT)-based method for eval- 
uating an F matrix. The F matrix is then supplemented 
with self-energy matrices £i ; 2 describing an open system 
connected to the two contacts and involving the detailed 
contact microstructure, while the noncquilibrium (trans- 
port) part is set up using the noncquilibrium Green's 
function (NEGF) formalism [27]. We refer the reader 
to our past work for details of the NEGF equations and 
the calculation of the self-energy matrices [28,29,17]. In 
this article, we will concentrate mainly on the physical 
insights. 



Poisson-hydrodynamic equations used extensively in ana- 
lyzing device transport [30,31]. We supplement the Pois- 
son (Hartree) term with exchange-correlation corrections 
that are small in macroscopic devices but significant in 
molecules, while the semiclassical hydrodynamic equa- 
tion is replaced with a fully quantum mechanical NEGF 
formalism [32]. 

The self-consistent formalism described above is com- 
pletely general, and can be employed in modeling trans- 
port through a wide variety of physically different sys- 
tems. For instance, we have used this scheme to obtain 
semi-empirical [3,33] and first-principles density func- 
tional (DFT)-based (LANL2DZ/B3PW91) [17,36] I-V 
characteristics of metallic quantum point contacts as well 
as semi-conducting aromatic thiol molecules bonded to 
Au(lll) surfaces. Whether one is dealing with a carbon 
nanotube, a ballistic MOSFET [33,34], a spin transis- 
tor, a resonant tunneling diode [35] or a molecular wire, 
the above formalism holds. One needs simply to evalu- 
ate each of the following quantities: (i) an appropriate 
Fock matrix F describing the device; (ii) a contact Fermi 
energy Ep relative to which the device levels are eval- 
uated and which determines the electrochemical poten- 
tials jt/1,2; (iii) a self-consistent potential J7scf describing 
charging effects and the electrostatic influence of the en- 
vironment (this term is incorporated into the effective F 
matrix of the device) , and (iv) a set of self-energy matri- 
ces Ei,2 that describe the coupling of the device with the 
contact, the matrices depending on the contact surface 
Green's function as well as the device-contact bonding 
geometry. Additional self-energy matrices can be intro- 
duced to describe scattering, by phonons or polarons for 
example. Within the same self-consistency scheme and 
NEGF prescription, one can then get qualitatively differ- 
ent I-V characteristics just by varying F, Ep, Uscf and 
S. 
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FIG. 2. Onset of conduction is given by the voltage where 
either of the contact chemical potentials /Lti,2 crosses the near- 
est conducting molecular level, HOMO in this case. 

The coupled DFT-NEGF formulation of molecular 
electronics is in effect the generalization of the coupled 
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FIG. 3. Self-consistency schemes: (a) for an isolated 
molecule in equilibrium, one calculates the Fock matrix F 
starting with a guess density matrix p, and fills up the cor- 
responding levels with N electrons to get back p; (b) for an 
open system, the molecular Fock matrix is supplemented with 
self-energy matrices £1,2 describing coupling with the con- 
tacts. An applied bias drives the system out of equilibrium 
due to two different contact chemical potentials /ii,2- The 
step from F to p is different from (a), and is obtained by 
solving the NEGF equations [32] . 
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The effect of each of the above parameters on the 
molecular I-V is schematically shown in Fig. 1. The con- 
ductance gap depends on (Ep — Smol), and the max- 
imum current level is set by the parallel combination 
r o ff = rir 2 /(Ti +r 2 ) of the individual contact broaden- 
ings I\ 2 = «(£i,2 — S| 2 ). The current rises over a width 
which depends on the total broadening of the molecular 
levels, which in turn depends on (a) the thermal broad- 
ening fcflT, (b) the series combination T = Ti + T 2 of the 
contact broadenings Ti j2 and (c) the Coulomb charging 
energy Uo to add an extra charge to the molecule (the 
charging energy tends to drag out the conductance peak, 
so that for appreciable charging energies, the contact 
chemical potentials may not be able to cross a molecular 
level easily under bias, resulting in a relatively featureless 
I-V characteristic). Finally, the potential profile across 
the molecule sets the overall voltage division factor r\ [3] 
which determines the prefactor in the ratio between the 
conductance gap and Ep — -Emol {v = 0.5 and the pref- 
actor equals 4 for symmetric coupling, as in Fig. 1). The 
voltage division factor r\ depends on the contact geome- 
tries and characterizes the Laplace part of the potential, 
while the Poisson part describes self-consistent charging 
effects, and is characterized by Uq. 

We summarize below the most challenging and phys- 
ically relevant questions in obtaining a molecular I-V 
characteristic: 

• Where is the contact Fermi energy relative 
to the molecular levels? (Ep, Emol) 

• What is the broadening due to the contacts? 

(r\ 2 ) 

• What is the spatial profile of the Laplace po- 
tential? (r/) 

• What is the charging energy? (Uq) 

Note that each of the above quantities is in general a 
complicated matrix that can be modeled independently 
using either semi-empirical or first-principles methods. 
However, in order to develop a "feel" for how these quan- 
tities affect the molecular I-V as in Fig. 1, we will try 
to capture their essence in terms of a few characteristic 
scalar parameters, as defined above. We will now address 
the influence of each parameter one at a time below. 



II. MOLECULAR CONDUCTION: HOW CAN WE 
MODEL IT? 

A. Where is the contact Fermi energy relative to the 
molecular levels? 

This is probably the most challenging problem to sort 
out. One needs to start by modeling the quantum chem- 
istry of the molecule. For our candidate molecule phenyl 



dithiol (PDT) , shown in Fig. 4, it is believed that the sul- 
phur atoms bond with a Au(lll) surface by desorption 
of the end hydrogen atoms that are then replaced by a 
triangle of gold atoms with sulphur sitting above their 
centroid [37]. The energy levels of the isolated PDT 
molecule compare well with those obtained by replac- 
ing each H-atom by three gold atoms (the gold atoms 
introduce some additional localized levels in the HLG). 
Replacing the gold cluster with a self-energy describing 
metallic gold broadens the molecular levels into a quasi- 
continuous spectrum, with the localized levels developing 
into metal-induced gap states (MIGS) decaying spatially 
away from the contacts into the molecular center. Now, 
the energy levels of the isolated molecule itself depend 
sensitively on the method of calculation (a comparison 
plot is shown in Fig. 4). Different theoretical groups have 
adopted different first-principles schemes in their analy- 
sis [17-20], so the unanswered question at this point is: 
which method is appropriate for calculating the single- 
particle energy levels of an open subsystem under bias ? 

PDT Energy levels in LANL2DZ basis set 
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FIG. 4. The single-particle molecular energy levels of 
PDT vary considerably for density functional-based and 
Hartree-Fock methods, giving different conductance gaps. 
The orbital wave functions agree for all the different meth- 
ods and for various choices of basis sets (symbols H: HOMO, 
L: LUMO, HLG: HOMO-LUMO gap, LDA: Local Density 
Approximation, B3PW91: 3 parameter Becke exchange and 
Perdew-Wang 91 correlation, HF: Hartree-Fock ). 

What is consistent among the various methods of cal- 
culation, including semi-empirical (Hiickel-based) theo- 
ries, is the overall shape of the orbital wave functions; 
for instance, the HOMO is largely sulphur-based and de- 
localized while the LUMO is ring-based and localized. 
This would seem to suggest a broad HOMO DOS and 
a sharp LUMO DOS, although theories don't seem to 
agree even on this point (see for example Fig. 22 in [28] 
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and Fig. 3 in [38]). The difference could arise due to 
different bonding geometries assumed at the gold-thiol 
interface. Thus the intrinsic chemistry of the molecular 
system needs to be cleared up, and theoretical agreement 
reached on the equilibrium molecular bonding properties 
before the transport-related issues can be sorted out. 

The position of the contact Fermi energy Ep relative 
to the molecular energy levels is also an unsettled is- 
sue. The Fermi energy depends sensitively on the spe- 
cific model for the contact geometry. Different models for 
the contact (for instance DFT-based [20,17], Bethe lat- 
tice [21], jcllium [19]) can give different charge contents 
and level broadenings of the molecule, as well as different 
work functions for the bulk electrodes. This could cause 
an appreciable shift in Ep, given the rather small gap 
density of states in the molecular HLG. Given that for 
extensively studied systems such as metal-semiconductor 
interfaces the precise location of Ep is still an active topic 
of research [39] , perhaps the best one can do at this point 
is to inquire if Ep is closer to the molecular HOMO or 
LUMO level, the analogous question for a semiconductor 
being whether it is p-type or n-type. 

Conceptually the cleanest way to address the equilib- 
rium Fermi energy problem is by including a few layers 
of the contact as a cluster in a "supermolecule" , which 
would act as the device under investigation. Thereafter 
the contact is assumed unaltered during conduction, with 
all the "action" lying in the device sector. The ad- 
vantages of including such a cluster are enormous (for 
a discussion, see [36]), such as the automatic inclusion 
of image charges (the supermolecule is charge neutral), 
avoiding uniqueness issues related to partitioning in a 
non-orthogonal, non-localized, atomic basis set and the 
proper treatment of the surface physics. Ideally, the clus- 
ter size should be significantly larger than the atomic 
Debye length of the contact material, while for practi- 
cal purposes, it is usually limited by computational re- 
sources. Ep is usually set by the HOMO of the (large) 
supermolecule, while the molecular levels can be identi- 
fied by either plotting the wavefunctions or by computing 
the local density of states (LDOS) on the molecule. To 
employ this scheme to sort out the molecular chemistry 
and energy level structure, it is essential that the contact 
cluster and the molecule be calculated using the same 
scheme (DFT/tight-binding, etc). Attempts at perform- 
ing such a computation at the semi-empirical [10] and 
DFT [40] levels have yielded a Fermi energy quite close 
to the HOMO level for PDT-Au(lll) hctcrostructures. 

The experimental situation is somewhat unclear. The 
conductance gap for PDT itself is different for differ- 
ent experimental geometries. The gap is around 3 volts 
for breakjunction measurements by Reed et al. [1] and 
around 4 volts for STM measurements by Hong et al. 
[41], while corresponding breakjunction measurements by 
David Janes' group at Purdue indicate featureless I-V 
characteristics with no discernible conductance gap [42] . 



Since the gap depends on Ep which is likely to be dif- 
ferent for the two experimental contact geometries, and 
could further be compromised by the presence of charg- 
ing in the system, such a difference is not surprising. It is 
also not clear whether the conduction is through a single 
molecule bridging the junctions, or a series combination 
of molecules attached separately to the two junctions [43] . 
It seems sensible therefore to treat Ep as a fitting pa- 
rameter, in the absence of precise characterization of the 
contact surfaces and molecular geometry. Alternatively, 
one could dictate the position of Ep relative to the lev- 
els, guided by separate equilibrium cluster calculations, 
as discussed earlier. 

Experiments incorporating a third terminal (gate) can 
help clarify some of these issues appreciably. For in- 
stance, a positive gate voltage lowers the molecular lev- 
els so that the Fermi energy approaches the LUMO and 
moves away from the HOMO. For a purely electrostatic 
gate control mechanism, a measured decrease in conduc- 
tance will suggest HOMO (p-type) conduction, while the 
reverse result indicates LUMO (n-type) conduction. 

Summary: Need to model the molecule and the contact 
bonding self- consistently within the same scheme, doing 
justice to the molecular quantum chemistry as well as 
the contact surface micro structure. The method of cal- 
culating the energy levels or the Fermi energy is still an 
unresolved issue. 

B. What is the broadening due to the contacts? 

Although experimental knowledge of the contact con- 
ditions is difficult to access, could one at least hope to 
model a particular idealized contact geometry and obtain 
an appropriate self-energy? We obtain the self-energy 
matrices Si. 2 formally by an exact partitioning of the in- 
finite metal-molecule-metal system, projecting its single 
particle Green's function onto the device subspace [36]. 
Si ; 2 depend on the contact surface Green's function and 
the contact-molecular bondings. We obtain the couplings 
at the surface by simulating a large cluster from the con- 
tact coupled to the device and calculating its overlap and 
Fock matrices. The contact surface Green's function is 
calculated using a recursive Green's function technique 
taking the full group theory of the FCC Au(lll) crystal 
surface into account [11,17]. 

One can replace the partitioning scheme with a 
scattering-formalism [19,12,10,13] that deals with the en- 
tire infinite system. Ideally, both methods (scattering 
formalism and the NEGF prescription) should yield the 
same answer; however, there is a conceptual simplifica- 
tion in partitioning the problem into a "device" part in- 
volving the electronically active molecule, and a "con- 
tact" part determining the lead-molecular interactions. 
These involve two entirely different areas of research, 
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quantum chemistry and surface physics, so partitioning 
allows us to improve modeling each of them indepen- 
dently. Moreover, NEGF naturally allows us to include 
incoherent processes, which can be important even for a 
short molecule if there are localized states that cannot 
be populated from the contacts [28,36]. 

How do we know if we have modeled our contacts cor- 
rectly? The overall "shape" of the molecular I-V can 
be obtained approximately without getting the bonding 
or the quantum chemistry right. One excellent bench- 
mark is the quantum point contact (QPC) [23], the I- 
V characteristic of which is experimentally measured 
to be ohmic, with a conductance quantized in units of 
Go = 2e 2 /h w 77 'fiS [44]. Starting with a six atom gold 
chain coupled to Au(lll) contacts describing a QPC, 
we get a conductance quantized ohmic I-V, as expected 
[17,36]. This is highly nontrivial, because conductance 
quantization arises out of a molecule that is perfectly 
transmitting over a band of energies between /ii and /12 
(aside from Fabry-Perot type oscillations). This requires 
the self-energy matrices to couple the wire with the con- 
tacts seamlessly without introducing any spurious reflec- 
tions. To illustrate the sensitivity of the quantization on 
the coupling, we scaled the overall matrix elements by a 
factor of five; the resultant I-V ceases to be ohmic, and 
resembles that of a resonantly conducting system such as 
PDT (Fig. 5 (i)). 
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FIG. 5. Calculated two-terminal I-V characteristics for dif- 
ferent molecular geometries: (i) ohmic I-V with a quantized 
conductance (adapted from [17]) for a quantum point con- 
tact (QPC) consisting of a six-atom gold chain connected to 
Au(lll) contacts; (ii) a symmetric, resonant I-V for PDT that 
turns asymmetric (upper and lower curves) (adapted from 
[57]) on altering the relative coupling strengths to the con- 
tacts; (iii) negative differential resistance (NDR) in the I-V 
for a QPC with a barrier in the middle (adapted from [36]). 

The above exercise is a good check of the accuracy 
of the contact broadenings. Once the surface Green's 
functions are deemed to be correct, it is an easy matter 
to replace the Aug molecular cluster with the molecule 
of choice and proceed with calculating its I-V. 



Summary: The QPC can be used as a benchmark for 
testing out the self-energy matrices. The couplings at the 
surface and the contact surface Green's functions need to 
be calculated accurately, including the overall group the- 
ory of the metal crystal comprising the leads. 

C. What is the spatial profile of the Laplace 
potential? 

The electrostatic potential profile across the molecule 
can be separated into two parts: (i) the Laplace part 
describes the influence of the contact geometries in the 
absence of charges and screening effects; (ii) the Poisson 
part involves screening by the charges, and is determined 
by the charging energy of the molecule. In this section, 
we will concentrate on the Laplace part, and address the 
charging-related issues associated with the Poisson part 
in the next section. 

The Laplace part of the potential profile is set by the 
relative capacitances of the contacts, and can be de- 
termined rigorously by solving the 3-D Laplace equa- 
tion with the correct potential boundary conditions for 
the contacts. The Laplace part is characterized by the 
voltage-division factor rj [3] which describes the propor- 
tion in which the applied voltage drops across the var- 
ious contact-molecular interfaces. A convenient way to 
analyze the potential profile across the molecule is to in- 
corporate the voltage-division factors in the definitions of 
the contact electrochemical potentials /ii,2, which would 
keep the molecular levels themselves fixed under drain 
bias in the absence of charging (next section) and move 
them under gate bias alone. For a two-terminal device 
appreciable current flow requires the capacitive couplings 
with the source and drain electrodes to be roughly equal 
(although their resistive (quantum) couplings Ti.2 could 
still be quite different), leading usually to rj « 0.5. 

Although the Laplace part involves essentially nine- 
teenth century classical electrostatics, it can substantially 
influence device I-V characteristics. In a gated three- 
terminal device, for instance, a good gate control mecha- 
nism in a well-designed ballistic conductor essentially in- 
volves trying to keep the charge density near the source 
end of the conductor constant by pinning the device DOS 
to the source chemical potential [51]. This gives an ef- 
fective 77 <C 1, so that as the drain chemical potential 
ventures into the HLG under bias, the current starts to 
saturate [52] . Note that this saturation mechanism is en- 
tirely different from saturation in two-terminal molecular 
devices, which occurs when either contact chemical po- 
tential has just crossed a molecular level and another level 
has not yet kicked in. In contrast to the two-terminal I- 
V, the three-terminal characteristic is asymmetric with 
respect to source-drain bias. Since the gate determines 
the position of the equilibrium Fermi energy through the 
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Laplace solution, it leads to gate-modulation of the ON 
current in such devices. 

A gate can influence the molecular electronic proper- 
ties in a variety of ways: it could have a purely electro- 
static effect on the channel charge and correspondingly 
the molecular levels, as described above. However, as 
pointed out by Damle et al. [52], good electrostatic gate 
control is not possible for a 10 A molecule such as PDT, 
unless the gate oxide is prohibitively thin. Additionally, 
the gate can alter the properties of the contact-molecular 
interfaces (Schottky-barrier type effects [53]), or even al- 
ter the conformations of the active molecule [54], all of 
which could affect the shape of the molecular I-V. 

Another example where the Laplace solution it- 
self can influence the molecular I-V characteristic in- 
volves conduction mechanisms that require the align- 
ment/misalignment of energy levels localized on different 
parts of the molecule. For instance, the Aviram-Ratncr 
mechanism [49] involves a donor-bridge-acceptor system 
where the Laplace part of the self-consistent potential 
aligns the levels at the two ends for positive bias, and 
misaligns them for negative bias, leading to a strongly 
rectifying I-V characteristic. A similar example involves 
a quantum point contact (QPC) with a stretched bond 
in the middle (Fig. 5(iii)). The defect disconnects the 
LDOS on its two sides, allowing them to separately equi- 
librate with the two contacts and follow their respective 
chemical potentials under source-drain bias. Within the 
window set by \i\ and /X2, the Laplace solution causes 
the LDOS on the two sides to slide past each other (fig- 
ure inset). Since this amounts to two transmission peaks 
sliding in and out of resonance, the resulting I-V shows 
a weak negative differential resistance [36,50]. 
Summary: The Laplace part of the 3-D electrostatic po- 
tential profile needs to be calculated using the boundary 
conditions set by the electrodes. The Laplace potential 
can significantly affect the molecular I-V characteristic, 
by aligning or misaligning different parts of the molecu- 
lar LDOS or by pinning the molecular DOS to one of the 
electrodes in the presence of a gate terminal. 

D. What is the charging energy? 

The Poisson part of the potential profile carries infor- 
mation about screening and charging inside the molecule, 
and is characterized by UqDq — Cq/Cs, where Uq = 
e 2 /Cs is the Coulomb charging energy, Do is the molec- 
ular gap DOS, the quantum capacitance Cq = e 2 D , 
and Cs describes the total electrostatic capacitive cou- 
pling to the various electrodes. The charging energy 
Uq describes the "ease" with which the molecular lev- 
els can be filled or emptied, and tends to drag out the 
I-V characteristic as shown in Fig. 1. Note that the effect 
of charging can be looked upon as a voltage-dependent 



r)(V). The net capacitance Cs contributing to Uq is de- 
termined by the geometry and dielectric constant of the 
molecule and the electrodes. For instance, while 10 A 
InAs spherical quantum dots have low charging energies 
<~ 100 meV [55], a 10 A isolated PDT molecular wire has 
a much larger charging energy ^3 — 4 cV, making it a lot 
harder to cross a level with a drain bias. At even higher 
charging energies Uo Li,2, one can get many-body ef- 
fects such as Coulomb Blockade and Kondo resonance. 
While some Coulomb-Blockade type effects can be cap- 
tured within an effective one-particle self-consistent field 
model within an unrestricted calculation [29] (such as 
unrestricted Hartree-Fock or spin density functional the- 
ory), doing justice to these problems requires us to go 
beyond the one-particle prescription we have used so far. 

The Poisson solution describes the effect of adding or 
removing charge from the molecule, as well as the effect 
of redistribution of charge within the molecule, respon- 
sible for screening of the applied voltage. The efficiency 
of the screening process depends on the amount of ma- 
terial available for the reorganization of the charges. For 
a molecular wire much thinner than the Debye length, 
the net electrostatic potential is essentially given by the 
Laplace solution of the previous section, leading to a 
ramp-like potential profile [17,36,45-47]. In contrast, a 
thick metallic wire allows sufficient screening, yielding a 
potential profile that is essentially flat [36,45,46]. Such 
a flat potential profile can be obtained even with a thin 
molecular wire if the latter is embedded in a dielectric 
medium as part of a self-assembled monolayer (SAM), 
in which case the neighboring wires screen the potential 
profile. One way to model a SAM would be disallow any 
transverse variations in charge or electrostatic potential, 
which would amount in effect to solving the 1-D Poisson 
equation. This would give a highly screened potential 
profile even with a thin molecular wire [48,45]. 

Charging can lead to very interesting effects, such as 
the creation of an asymmetric I-V with a spatially sym- 
metric molecule [5,56]. Consider a symmetric molecule 
with unequal resistive (quantum) couplings to the two 
contacts (Ti ^ I^). Near the onset of current conduc- 
tion through a HOMO level, a negative bias on the strong 
contact keeps it filled, while a positive bias empties it. 
Since the molecule gets positively charged one way but 
not the other, the I-V is dragged out asymmetrically such 
that one gets a lower current for positive bias on the 
stronger contact. Interestingly for conduction through 
a LUMO level, the sense of the I-V asymmetry reverses 
[57], because one now needs to fill the LUMO to charge 
up the molecule. This allows us to identify the nature 
of the conducting molecular orbital, which is important 
given that different orbitals conduct quite differently. For 
PDT, STM data [3] (the STM tip being the weaker con- 
tact) seems to indicate HOMO-based conduction. Our 
results (Fig. 5 (ii)) qualitatively match the I-V charac- 
teristics obtained by Reichert et al. [5], with an initially 
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symmetric I-V that turns weakly asymmetric on drawing 
either contact away from the molecule. 

The total electrostatic potential profile is the sum 
of the Laplace and Poisson parts, which have differ- 
ent effects on the molecular I-V. For typical break- 
junction/STM measurements on molecules such as PDT, 
Cq/Cy, = UqDo <C 1, and the overall potential profile 
is essentially given by the Laplace potential with the net 
voltage dividing up as the capacitance ratio r\ between the 
electrodes. If however the molecule is so strongly coupled 
to the substrate that there are appreciable MIGS, D be- 
comes large enough that Cq/Cs starts to become impor- 
tant. In the limit of Cq/Cy, 3> L the I-V gets dragged 
out substantially by charging (Fig. 1), and the final self- 
consistent voltage divides between the source and drain 
contacts according to their resistance ratio T^/ (Ti + T%). 
This situation arises for metal nanoclusters probed by an 
STM tip. Since the resistance is much larger at the STM 
end, much of the applied voltage drops across the STM- 
cluster gap, so that the cluster levels remain pinned to the 
substrate, and the STM chemical potential alternately 
scans the HOMO and LUMO levels under opposite drain 
bias. 

Summary: The charging energy can turn an otherwise 
symmetric I-V into an asymmetric one. Given a spa- 
tially symmetric molecule, we predict a larger current for 
positive bias on the stronger contact and HOMO conduc- 
tion, while for LUMO conduction, the sense of the I-V 
asymmetry is reversed. 

III. COMPUTATIONAL ISSUES 

A brief discussion of computational issues is perhaps in 
order. The calculation of molecular conductance requires 
two steps, (a) calculating a Fock matrix given a density 
matrix, and (b) calculating a density matrix from the 
Fock matrix. The first step is the most computationally 
challenging part, involving the evaluation of DFT-based 
exchange-correlation integrals which are quite numeri- 
cally complicated, especially in sophisticated basis sets 
involving rclativistic core pseudopotentials. We find our 
own LDA calculation of the molecular Fock matrix in a 
minimal basis set to be comparable in accuracy, but not 
in speed, with GAUSSIAN98. Therefore we let GAUS- 
SIAN98 do this part, exploiting decades of development 
that have gone into it. The second step requires contact 
self-energy matrices which can be calculated once and for 
all for realistic contact surface structures (e.g., Au(lll) 
contacts) using a real space recursive Green's function 
technique. The contact coupling matrices can be simu- 
lated in GAUSSIAN 98 with a finite-sized cluster, with 
care exercised to eliminate edge effects on the structure 
of the contact surface Green's function [36] (a suitable lo- 
calized basis describing gold would sort out this problem 



automatically). 

The computational challenge for us is to solve the 
NEGF equations, requiring us to find the number of elec- 
trons on the molecule. Such a requirement amounts to 
integrating the nonequilibrium electron DOS all the way 
from the bottom of the contact band to the Fermi en- 
ergy. Since the Green's functions entering the DOS ex- 
pression are highly peaked around the molecular levels, 
this process involves integrating over an energy range of 
several hundred volts with a millivolt accuracy for each 
bias point and each step of the self-consistent procedure. 
We have addressed this problem in two ways: (i) assum- 
ing a weak (in practice, constant) energy-dependence of 
the self-energy matrices (valid for Au(lll) which has an 
energy-independent DOS near Ep), one can then per- 
form the integrals analytically [36]; (ii) the nonequilib- 
rium DOS can be divided into two groups, one lying out- 
side the domain of //i. 2 , and the rest inside that window. 
The first part can be integrated using a contour integra- 
tion scheme [20], while the part between fii^ can be cal- 
culated either by brute force grid-based integration over 
the finite range between ji\ and /i 2 , or by reverting to the 
constant S approximation. Although our calculations are 
performed with LANL2DZ basis sets which are somewhat 
delocalized, it is preferable to employ relatively localized 
basis sets in order to avoid issues related to partition- 
ing and to be consistent with the tight-binding approach 
that we are using here [23] . Different approximations give 
different values of the total electron count, thereby affect- 
ing the equilibrium Fermi energy position, which clearly 
needs more attention. However, the approximations al- 
low us to obtain first-principles DFT-based I-Vs for a 
molecule like PDT in a few hours on a SUN workstation, 
taking into account the details of the contact geometry. 

IV. ISSUES WE HAVEN'T COVERED 

We have seen that by appropriately modeling each ex- 
perimental geometry, we can get qualitatively and quan- 
titatively different conductance properties, ranging from 
ohmic to rectifying, switching and saturating I-V char- 
acteristics. We now outline three of the issues that we 
have ignored so far, namely, Conformation, Incoherence 
and Correlation. 

(i) Conformation. One of the principal advantages of 
a molecule is that it is semi-flexible. This means that the 
conformation of a molecule can be altered, by transfer- 
ring charge or applying an external field. Encouraging 
experimental indications of a conformationally mediated 
I-V have been obtained, for a fullerene-based transistor 
[58], and for the redox sidegroup-specific [2] or vibra- 
tionally mediated [59] NDR measurements. We are cur- 
rently investigating the role of conformational changes, 
in conjunction with charging and gate electrostatics, in 
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modulating the molecular I-V [54] in a three-terminal de- 
vice. 

(ii) Incoherence. For long molecules, molecular vibra- 
tions are important as sources of incoherent or inelastic 
scattering, leading to hopping or polaronic transport (see 
for e.g. [60]). Such inelastic effects can be naturally in- 
cluded in the NEGF formalism, requiring us to introduce 
another self-energy matrix describing the connection of 
the molecule with the source of the inelastic scattering 
(a phonon bath, for example) [36,45,61]. In the tunnel- 
ing regime, in particular, inelastic scattering turns out to 
be crucial for transport and dissipation in long molecular 
wires, such as DNA. 

(iii) Correlation. Finally, there are examples of molec- 
ular measurements such as the Kondo effect [62], where 
the physics of current flow cannot be captured in terms 
of a simple one-particle picture, and require the incorpo- 
ration of many-body correlation effects in our model. We 
leave these problems for future work. 

We would like to thank P. Damle, D. Janes, A. Liang, 
S. Lodha, M. Lundstrom, M. Paulsson, T. Rakshit, R. 
Reifenberger, R. Venugopal and F. Zahid for useful dis- 
cussions. This work has been supported by the US Army 
Research Office (ARO) under grant number DAAD19- 
99-1-0198. 
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